肉质性状表型矫正

#garden/🌱Seedling

早期纠结的问题

需要解决的问题：

是否存在、不同实验室/仪器效应、不同品种、不同个体（年龄、性别）
- 四个实验室都测的几个品种，同一品种内部，检测是否有实验室差异
- ~~所有样本，看实验室是否差异（看其他因素是否可以作为固定效应）~~不太可行，品种差异
- 最多的两家单位内部、同一品种内部，是否存在批次效应
- 其实年龄和性别也会影响物理性质，但先假设均匀分布，不去考虑，后期专门拿出来看。
如何矫正不同实验室不同仪器的系统误差
- 标准样品才可以进行
- 拟合线性模型，分离不同实验室造成的系统误差
[ ]

下一步：

月龄260+？
批次效应
矫正后的定性和定量分析是否准确性高一些
可以先GWAS

矫正值 = 固定效应预测值 + 残差

啥意思，残差不就是矫正后的表型值吗，为什么要固定效应加上残差？

表型矫正之后，为什么数据分布更集中了

关于表型矫正的思考

什么时候需要表型矫正？

Question

还是有问题想请教你，我很疑惑光谱数据拟合模型的时候要不要表型矫正。

比如我们GS预测出来的其实不是实际的表型值，而是GEBV，通过基因组数据只能预测基因组决定的这部分，真实的表型值没法预测。

同理，我们通过光谱想预测肉类的营养物质含量指标或者奶牛乳成分之类的，参考个体是用矫正后的表型值去训练，然后拿光谱数据一起构建模型，然后通过光谱去预测未知个体的营养物质之类的表型，其实也无法预测实际表型吧，预测出来的也是类似矫正后的表型值是吗，其实还是不知道这个个体真实的营养物质含量。这样的话，是不是在实际生产当中就意义不大了，因为其它因素没有办法表征，预测值和真实值之间有相关性但总会有差异，只能想GEBV一样看个排名？

还是说用光谱预测直接不需要矫正，就拿实际表型值训练，然后预测出来也是实际表型值，但我觉得不矫正准确度不会太高啊，因为总会有其它因素去干扰。你们之前光谱数据表型值需要矫正吗

Tldr

做GS和GWAS需要表型矫正，做光谱数据定量模型的时候不需要表型矫正。

Success

因为GWAS和GS就是要去除其它效应，只把遗传效应纳入考虑，衡量遗传效应。
而通过光谱预测的时候，就是想直接得到物质的浓度或者含量，如果使用矫正表型进行拟合，那么就无法预测实际的含量了。

矫正表型应该用残差，均值加残差，还是固定效应预测值加残差

矫正值 = 固定效应预测值 + 残差

啥意思，残差不就是矫正后的表型值吗，为什么要固定效应加上残差？

Success

关键还是要看用途，如果是用来GWAS，那直接用残差就可以，或者整体均值加残差，都一样的。
如果是想保留品种间的差异，但我们矫正的时候把品种作为固定效应了，这时候就需要固定效应的的预测值加残差，相当于再把品种这部分效应再加上，这样整体的数据才会有品种间的差异，如果不加的话，相当于直接矫正掉了。

表型矫正代码python

R语言中lm和glm与lmer的区别

statsmodels库中sm.OLS和smf.ols的区别

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

def adjustment(df, num_covs_lst, class_covs_lst, pheno):
    """
    表型矫正函数
    
    参数:
    df: 包含所有数据的DataFrame
    num_covs_lst: 数值型协变量列表 (如 ['月龄'])
    class_covs_lst: 分类协变量列表 (如 ['性别', '品种名称'])
    pheno: 要矫正的表型列名 (字符串)
    
    返回:
    result: 模型拟合结果
    residuals: 残差（矫正后的表型值）
    """
    # 1. 构建公式字符串
    formula_parts = []
    
    # 添加数值型协变量
    if num_covs_lst:
        num_part = " + ".join(num_covs_lst)
        formula_parts.append(num_part)
    
    # 添加分类协变量
    if class_covs_lst:
        class_part = " + ".join([f"C({var})" for var in class_covs_lst])
        formula_parts.append(class_part)
    
    # 如果没有协变量，只使用截距
    if not formula_parts:
        formula_parts = ["1"]
    
    # 组合完整公式
    formula = f"Q('{pheno}') ~ " + " + ".join(formula_parts)
    
    # 2. 创建模型并拟合
    try:
        model = smf.ols(formula, data=df)
        result = model.fit()
        
        # 3. 提取残差
        residuals = result.resid
        
        # 4. 创建完整的残差序列（包含原始索引）
        full_residuals = pd.Series(np.nan, index=df.index)
        valid_index = result.resid.index  # 模型实际使用的索引
        full_residuals.loc[valid_index] = residuals
        
        return result, full_residuals
    
    except Exception as e:
        print(f"Adjustment failed: {str(e)}")
        print(f"Formula: {formula}")
        return None, pd.Series(np.nan, index=df.index)

使用方法：

# 矫正表型
df=pd.read_csv('meat_quality_factors.csv')
phenotype_columns = ['肉色L', 'pH值', '水分']
# 定义协变量
num_covs = ['月龄']  # 数值型
class_covs = ['性别', '品种名称']  # 分类型
# 批量矫正
corrected_data = df.copy()
for pheno in phenotype_columns:
    _, residuals = adjustment(df, num_covs, class_covs, pheno)
    corrected_data[f'adjusted_{pheno}'] = residuals



# 保存矫正后的表型作为GWAS输入
#gwas_input = corrected_data[['样品编号'] + [f'adjusted_{pheno}' for pheno in phenotype_columns]]
#gwas_input.columns = ['样品编号'] + phenotype_columns
#gwas_input
# 保存到文件
#gwas_input.to_csv('adjusted_phenotypes.csv', index=False)

# 保存矫正后的表型作为GWAS输入
gwas_input = corrected_data[['实验室编号']].copy()  # 创建副本避免修改原始数据

# 添加每个矫正后的表型列，并显式设置NA值
for pheno in phenotype_columns:
    col_name = f'adjusted_{pheno}'
    # 将np.nan转换为'NA'
    gwas_input[col_name] = corrected_data[col_name].apply(lambda x: 'NA' if pd.isna(x) else x)

# 保存到文件（显式设置NA表示）
gwas_input.to_csv('adjusted_phenotypes2.csv', index=False, na_rep='NA')